Various number of packages which are required for the 3D map generation. These are installed and loaded into the program. These packages avail the functions such as Plotting of points, shading and adjusting the colors, conversion of GIS data, validating GDAL operations,etc.

library(rayshader)
library(magick)
library(sp)
library(raster)
library(scales)
library(rgdal)
library(rgeos)
library(av)

The GIS data for elevation is downloaded from Derek Watkin’s -“30 meter SRTM tile downloader”. The downloaded SRTM hgt data is converted to a matrix. The result of the plot is an elevation intensity profile of the tile.

elevation1 = raster::raster("D:/Assam maps/New Folder/everest/N27E086.hgt")
elevation2 = raster::raster("D:/Assam maps/New Folder/everest/N27E087.hgt")
everest_elevation = raster::merge(elevation1,elevation2)
height_shade(raster_to_matrix(everest_elevation)) %>%
plot_map()
[1] "Dimensions of matrix are: 7201x3601."

The coordinate reference of the elevation data matched to the imagery data


crs(everest_elevation)
CRS arguments: +proj=longlat +datum=WGS84 +no_defs 

The desired area whose map is to be plot is cropped in a rectangular shape by defining the co-ordinates of the bottom left and top right, diagonally opposite points. The longitude and latitude of the two points are found using “Google Maps” and then defined.

bot_lefty=86.856740
bot_leftx=27.69857
top_righty=87.308946
top_rightx=27.991321
distl2l= top_righty-bot_lefty  
scl=1-(20/(distl2l*111))
bottom_left = c(bot_lefty, bot_leftx)
top_right   = c(top_righty, top_rightx)
extent_latlong = SpatialPoints(rbind(bottom_left, top_right), proj4string=sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
e = extent(extent_latlong)
e
class      : Extent 
xmin       : 86.85674 
xmax       : 87.30895 
ymin       : 27.69857 
ymax       : 27.99132 
scl
[1] 0.6015529

Then, the elevation data is converted into Base R matrix

elevation_cropped = crop(everest_elevation, e)

everestel_matrix = raster_to_matrix(elevation_cropped)
[1] "Dimensions of matrix are: 1628x1054."

The geographical plot with 3d visualization

everestel_matrix %>%
 
  sphere_shade(texture = "bw") %>%
  plot_3d( everestel_matrix, windowsize = c(1100,900), zscale = 15, shadowdepth = -15,
           zoom=0.5, phi=45,theta=-45,fov=70, background = "#F2E1D0", shadowcolor = "#523E2B")
render_scalebar(limits=c(0,10,20),label_unit = "km",position = "S", y=50,scale_length = c(scl,1))
render_compass(position = "W" )
render_snapshot(title_text = "Mount Everest_geographical|Imagery: Landsat 8 | DEM: 30m SRTM",title_bar_color = "#1f5214", title_color = "white", title_bar_alpha = 1)

optional: Conversion of the map to a video footage of its planar rotation

angles= seq(0,360,length.out = 1441)[-1]
for(i in 1:1440) {
render_camera(theta=-45+angles[i])
  render_snapshot(filename = sprintf("evr_geo%i.png", i), 
                  title_text = "Everest_geographical| DEM: 30m SRTM",
                  title_bar_color = "#1f5214", title_color = "white", title_bar_alpha = 1)
}
rgl::rgl.close()
system("ffmpeg -framerate 60 -i evr_geo%d.png -vcodec libx264 -an Mount_Everest_geographical_1.mp4 ")

Himalayas_geographical

DQotLS0NCnRpdGxlOiAiRXZlcmVzdCAzRCBNYXAgR2VuZXJhdGlvbiINCm91dHB1dDoNCiAgaHRtbF9ub3RlYm9vazogZGVmYXVsdA0KICBodG1sX2RvY3VtZW50Og0KICAgIGRmX3ByaW50OiBwYWdlZA0KLS0tDQoNCg0KVmFyaW91cyBudW1iZXIgb2YgcGFja2FnZXMgd2hpY2ggYXJlIHJlcXVpcmVkIGZvciB0aGUgM0QgbWFwIGdlbmVyYXRpb24uIFRoZXNlIGFyZSBpbnN0YWxsZWQgYW5kIGxvYWRlZCBpbnRvIHRoZSBwcm9ncmFtLiBUaGVzZSBwYWNrYWdlcyBhdmFpbCB0aGUgZnVuY3Rpb25zIHN1Y2ggYXMgUGxvdHRpbmcgb2YgcG9pbnRzLCBzaGFkaW5nIGFuZCBhZGp1c3RpbmcgdGhlIGNvbG9ycywgY29udmVyc2lvbiBvZiBHSVMgZGF0YSwgdmFsaWRhdGluZyBHREFMIG9wZXJhdGlvbnMsZXRjLiANCmBgYHtyfQ0KbGlicmFyeShyYXlzaGFkZXIpDQpsaWJyYXJ5KG1hZ2ljaykNCmxpYnJhcnkoc3ApDQpsaWJyYXJ5KHJhc3RlcikNCmxpYnJhcnkoc2NhbGVzKQ0KbGlicmFyeShyZ2RhbCkNCmxpYnJhcnkocmdlb3MpDQpsaWJyYXJ5KGF2KQ0KYGBgDQpUaGUgR0lTIGRhdGEgZm9yIGVsZXZhdGlvbiBpcyBkb3dubG9hZGVkIGZyb20gRGVyZWsgV2F0a2luJ3MgLSIzMCBtZXRlciBTUlRNIHRpbGUgZG93bmxvYWRlciIuIFRoZSBkb3dubG9hZGVkIFNSVE0gaGd0IGRhdGEgaXMgY29udmVydGVkIHRvIGEgbWF0cml4LiBUaGUgcmVzdWx0IG9mIHRoZSBwbG90IGlzIGFuIGVsZXZhdGlvbiBpbnRlbnNpdHkgcHJvZmlsZSBvZiB0aGUgdGlsZS4NCmBgYHtyIGZpZzEsIGZpZy5oZWlnaHQgPSAxNSwgZmlnLndpZHRoID0gMTAsIGFsaWduPSAiY2VudGVyIn0NCmVsZXZhdGlvbjEgPSByYXN0ZXI6OnJhc3RlcigiRDovQXNzYW0gbWFwcy9OZXcgRm9sZGVyL2V2ZXJlc3QvTjI3RTA4Ni5oZ3QiKQ0KZWxldmF0aW9uMiA9IHJhc3Rlcjo6cmFzdGVyKCJEOi9Bc3NhbSBtYXBzL05ldyBGb2xkZXIvZXZlcmVzdC9OMjdFMDg3LmhndCIpDQpldmVyZXN0X2VsZXZhdGlvbiA9IHJhc3Rlcjo6bWVyZ2UoZWxldmF0aW9uMSxlbGV2YXRpb24yKQ0KaGVpZ2h0X3NoYWRlKHJhc3Rlcl90b19tYXRyaXgoZXZlcmVzdF9lbGV2YXRpb24pKSAlPiUNCnBsb3RfbWFwKCkNCg0KYGBgDQpUaGUgY29vcmRpbmF0ZSByZWZlcmVuY2Ugb2YgdGhlIGVsZXZhdGlvbiBkYXRhIG1hdGNoZWQgdG8gdGhlIGltYWdlcnkgZGF0YQ0KYGBge3J9DQoNCmNycyhldmVyZXN0X2VsZXZhdGlvbikNCmBgYA0KDQpUaGUgZGVzaXJlZCBhcmVhIHdob3NlIG1hcCBpcyB0byBiZSBwbG90IGlzIGNyb3BwZWQgaW4gYSByZWN0YW5ndWxhciBzaGFwZSBieSBkZWZpbmluZyB0aGUgY28tb3JkaW5hdGVzIG9mIHRoZSBib3R0b20gbGVmdCBhbmQgdG9wIHJpZ2h0LCBkaWFnb25hbGx5IG9wcG9zaXRlIHBvaW50cy4gVGhlIGxvbmdpdHVkZSBhbmQgbGF0aXR1ZGUgb2YgdGhlIHR3byBwb2ludHMgYXJlIGZvdW5kIHVzaW5nICJHb29nbGUgTWFwcyIgYW5kIHRoZW4gZGVmaW5lZC4NCmBgYHtyfQ0KYm90X2xlZnR5PTg2Ljg1Njc0MA0KYm90X2xlZnR4PTI3LjY5ODU3DQp0b3BfcmlnaHR5PTg3LjMwODk0Ng0KdG9wX3JpZ2h0eD0yNy45OTEzMjENCmRpc3RsMmw9IHRvcF9yaWdodHktYm90X2xlZnR5ICANCnNjbD0xLSgyMC8oZGlzdGwybCoxMTEpKQ0KYm90dG9tX2xlZnQgPSBjKGJvdF9sZWZ0eSwgYm90X2xlZnR4KQ0KdG9wX3JpZ2h0ICAgPSBjKHRvcF9yaWdodHksIHRvcF9yaWdodHgpDQpleHRlbnRfbGF0bG9uZyA9IFNwYXRpYWxQb2ludHMocmJpbmQoYm90dG9tX2xlZnQsIHRvcF9yaWdodCksIHByb2o0c3RyaW5nPXNwOjpDUlMoIitwcm9qPWxvbmdsYXQgK2VsbHBzPVdHUzg0ICtkYXR1bT1XR1M4NCIpKQ0KZSA9IGV4dGVudChleHRlbnRfbGF0bG9uZykNCmUNCnNjbA0KYGBgDQpUaGVuLCB0aGUgZWxldmF0aW9uIGRhdGEgaXMgY29udmVydGVkIGludG8gQmFzZSBSIG1hdHJpeA0KYGBge3IgZmlnNCwgZmlnLmhlaWdodCA9IDE1LCBmaWcud2lkdGggPSAxMCwgYWxpZ249ICJjZW50ZXIifQ0KZWxldmF0aW9uX2Nyb3BwZWQgPSBjcm9wKGV2ZXJlc3RfZWxldmF0aW9uLCBlKQ0KDQpldmVyZXN0ZWxfbWF0cml4ID0gcmFzdGVyX3RvX21hdHJpeChlbGV2YXRpb25fY3JvcHBlZCkNCg0KYGBgDQoNClRoZSBnZW9ncmFwaGljYWwgcGxvdCB3aXRoIDNkIHZpc3VhbGl6YXRpb24NCmBgYHtyIGZpZzcsIGZpZy5oZWlnaHQgPSAxMiwgZmlnLndpZHRoID0gOSwgYWxpZ249ICJjZW50ZXIifQ0KZXZlcmVzdGVsX21hdHJpeCAlPiUNCiANCiAgc3BoZXJlX3NoYWRlKHRleHR1cmUgPSAiYnciKSAlPiUNCiAgcGxvdF8zZCggZXZlcmVzdGVsX21hdHJpeCwgd2luZG93c2l6ZSA9IGMoMTEwMCw5MDApLCB6c2NhbGUgPSAxNSwgc2hhZG93ZGVwdGggPSAtMTUsDQogICAgICAgICAgIHpvb209MC41LCBwaGk9NDUsdGhldGE9LTQ1LGZvdj03MCwgYmFja2dyb3VuZCA9ICIjRjJFMUQwIiwgc2hhZG93Y29sb3IgPSAiIzUyM0UyQiIpDQpyZW5kZXJfc2NhbGViYXIobGltaXRzPWMoMCwxMCwyMCksbGFiZWxfdW5pdCA9ICJrbSIscG9zaXRpb24gPSAiUyIsIHk9NTAsc2NhbGVfbGVuZ3RoID0gYyhzY2wsMSkpDQpyZW5kZXJfY29tcGFzcyhwb3NpdGlvbiA9ICJXIiApDQpyZW5kZXJfc25hcHNob3QodGl0bGVfdGV4dCA9ICJNb3VudCBFdmVyZXN0X2dlb2dyYXBoaWNhbHxJbWFnZXJ5OiBMYW5kc2F0IDggfCBERU06IDMwbSBTUlRNIix0aXRsZV9iYXJfY29sb3IgPSAiIzFmNTIxNCIsIHRpdGxlX2NvbG9yID0gIndoaXRlIiwgdGl0bGVfYmFyX2FscGhhID0gMSkNCmBgYA0Kb3B0aW9uYWw6IENvbnZlcnNpb24gb2YgdGhlIG1hcCB0byBhIHZpZGVvIGZvb3RhZ2Ugb2YgaXRzIHBsYW5hciByb3RhdGlvbg0KYGBge3J9DQphbmdsZXM9IHNlcSgwLDM2MCxsZW5ndGgub3V0ID0gMTQ0MSlbLTFdDQpmb3IoaSBpbiAxOjE0NDApIHsNCnJlbmRlcl9jYW1lcmEodGhldGE9LTQ1K2FuZ2xlc1tpXSkNCiAgcmVuZGVyX3NuYXBzaG90KGZpbGVuYW1lID0gc3ByaW50ZigiZXZyX2dlbyVpLnBuZyIsIGkpLCANCiAgICAgICAgICAgICAgICAgIHRpdGxlX3RleHQgPSAiRXZlcmVzdF9nZW9ncmFwaGljYWx8IERFTTogMzBtIFNSVE0iLA0KICAgICAgICAgICAgICAgICAgdGl0bGVfYmFyX2NvbG9yID0gIiMxZjUyMTQiLCB0aXRsZV9jb2xvciA9ICJ3aGl0ZSIsIHRpdGxlX2Jhcl9hbHBoYSA9IDEpDQp9DQpyZ2w6OnJnbC5jbG9zZSgpDQpzeXN0ZW0oImZmbXBlZyAtZnJhbWVyYXRlIDYwIC1pIGV2cl9nZW8lZC5wbmcgLXZjb2RlYyBsaWJ4MjY0IC1hbiBNb3VudF9FdmVyZXN0X2dlb2dyYXBoaWNhbF8xLm1wNCAiKQ0KYGBgDQohW0hpbWFsYXlhc19nZW9ncmFwaGljYWxdKEM6L1VzZXJzL2FzdXMvRG9jdW1lbnRzL01vdW50X0V2ZXJlc3RfZ2VvZ3JhcGhpY2FsXzEubXA0KQ==